library(dplyr)
library(tidyr)
library(ggplot2)
library(viridis)
library(ggpubr)
library(cowplot)
library(gridExtra)
library(stringr)

if (!file.exists("output")) {
  dir.create("output")
}

# load data
# Note: Signature attributions can be found in the accompanying Supplementary Tables of the HNC manuscript

source("HNC_metadata_tidy.R")
source("data_handling.R")

sbsCOSMIC <- read.csv("Input_signature_attributions/Input_SBS_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")
dbsCOSMIC <- read.csv("Input_signature_attributions/Input_DBS_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")
idCOSMIC <- read.csv("Input_signature_attributions/Input_ID_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")

additional_df <- list(sbsCOSMIC, dbsCOSMIC, idCOSMIC)

cosmic <- additional_df[[1]]
for (n in 2:length(additional_df)) {
  cosmic <- cosmic %>% left_join(additional_df[[n]], by = "donor_id")
}

# merge with cosmic sigantures and edit factor labels
data <- data %>%
  left_join(cosmic, by = "donor_id") %>%
  mutate(tobacco = factor(tobacco, labels = c("Non-smoker", "Ex-smoker", "Current-smoker")))

# Select signatures present in >n% of samples
# Dichotomize the signatures (Presence (signatures>0)="1",absence (signatures==0)="0")
# If signature is present in more than 75% cases, dichotomize by the median
data <- dichotomizeSigs(metadata = data, sigsn = cosmic, sigcutoff = 2)

# Plot
for (sig in c("SBS4", "SBS92", "SBS_I", "DBS2", "DBS6", "ID3")) {
  plot_kw(data,
    signatures = sig, factors = "tobacco",
    pval_cutoff = 1, y_labs = "Mutation burden"
  )
  ggsave(paste0("output/Figure_3b_", sig, "_", Sys.Date(), ".pdf"),
    device = "pdf",
    width = 3.5, height = 4, dpi = 700
  )
}
